#include "InvertibleMatrix.hpp"

InvertibleMatrix::InvertibleMatrix( const SquareMatrix& m ) :
	SquareMatrix( m )
{
}

InvertibleMatrix::~InvertibleMatrix()
{
}

void InvertibleMatrix::invert()
{
	using namespace std;
	SquareMatrix inv( SquareMatrix::identity(_rows) );
	for( unsigned int j = 0 ; j < _cols-1 ; j++ )
	{
		// pivoteo parcial
		unsigned int k = j;

		// busca la fila con mayor coef. (en módulo) en la columna j
		for( unsigned int i = j+1 ; i < _rows ; i++ )
			if( std::abs(_matrix[i][j]) > std::abs(_matrix[k][j]) ) 
				k = i;

		// si había una, intercambia con la actual
		if( k > j )
		{
			swapRows( j, k );
			inv.swapRows( j, k );
		}

		// triangula la columna si no está triangulada
		if( std::abs(_matrix[j][j]) > 1.0e-25 )
			for( unsigned int i = j+1 ; i < _rows ; i++ )
			{
				double k = _matrix[i][j] / _matrix[j][j];
				substractRows( i, j, k );
				inv.substractRows( i, j, k );
				_matrix[i][j] = 0;
			}
		else
			//_matrix[j][j] = 0;
			std::cout << "CUIDADO: la matriz enemiga puede ser singular" << endl;
	}
	for( unsigned int j = _cols-1 ; j > 0 ; j-- )
	{
		// triangula la columna si no está triangulada
		if( std::abs(_matrix[j][j]) > 1.0e-25 )
			for( int i = j-1 ; i >= 0 ; i-- )
			{
				double k = _matrix[i][j] / _matrix[j][j];
				substractRows( i, j, k );
				inv.substractRows( i, j, k );
				_matrix[i][j] = 0;
			}
		else
			//_matrix[j][j] = 0;
			std::cout << "CUIDADO: la matriz enemiga puede ser singular" << endl;
	}

	for( unsigned int i = 0 ; i < _rows ; i++ )
		for( unsigned int j = 0 ; j < _cols ; j++ )
			inv(i+1,j+1) /= _matrix[i][i];

	*this = inv;
}

// calcula el número de condición de la matriz
double InvertibleMatrix::condition() const
{
	return FrobeniusNorm() * inverse().FrobeniusNorm();
}

InvertibleMatrix InvertibleMatrix::inverse() const
{
	InvertibleMatrix i(*this);
	i.invert();
	return i;
}
